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Abstract 1.5D Vlasov-Maxwell simulations arc employed to model electromag- 
netic emission generation in a fully self-consistent plasma kinetic model for the 
first time in the solar physics context. The simulations mimic the plasma emission 
mechanism and Larmor drift instability in a plasma thread that connects the Sun 
to Earth with the spatial scales compressed appropriately. The effects of spatial 
density gradients on the generation of electromagnetic radiation are investigated. 
It is shown that 1.5D inhomogeneous plasma with a uniform background mag- 
netic field directed transverse to the density gradient is aperiodically unstable 
to Larmor-drift instability. The latter results in a novel effect of generation of 
electromagnetic emission at plasma frequency. The generated perturbations con- 
sist of two parts: (i) non-escaping (trapped) Langmuir type oscillations which are 
localised in the regions of density inhomogencity, and arc highly filamentary, with 
the period of appearance of the filaments close to electron plasma frequency in 
the dense regions; and (ii) escaping electromagnetic radiation with phase speeds 
close to the speed of light. When density gradient is removed (i.e. when plasma 
becomes stable to Larmor-drift instability) and a low density super-thermal, 
hot beam is injected along the domain, in the direction perpendicular to the 
magnetic field, plasma emission mechanism generates non-escaping Langmuir 
type oscillations which in turn generate escaping electromagnetic radiation. It is 
found that in the spatial location where the beam is injected, the standing waves, 
oscillating at the plasma frequency, are excited. These can be used to interpret 
the horizontal strips (the narrowband line emission) observed in some dynamical 
spectra. Quasilincar theory predictions: (i) the electron free streaming and (ii) 
the beam long relaxation time, in accord with the analytic expressions, are 
corroborated via direct, fully-kinetic simulation. Finally, the interplay of Larmor- 
drift instability and plasma emission mechanism is studied by considering dense 
electron beam in the Larmor-drift unstable (inhomogeneous) plasma. The latter 
case enables one to study the deviations from the quasilinear theory. 
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1. Introduction 

The aim of this paper is to employ Vlasov-Maxwell simulations for the elec- 
tromagnetic wave generation by a super-thermal (0.2 — 0.5c), hot electron beam 
injected into the solar coronal magnetised plasma. Since such beams are thought 
to be responsible for the generation of type III solar radio bursts we start from a 
brief review of the previous relevant results. Although we stress that for the rea- 
sons described below, at this stage, the model presented here cannot be directly 
applied to the type III busts and it is perhaps more relevant for the interpretation 
of the narrowband line emission observations. The type III solar radio bursts are 
believed to come from the super-thermal beam of electrons that travel away 
from the Sun producing the observed electromagnetic (EM) radiation via the 
plasma emission mechanism (see e.g. IMelrosel INindos e t al.\ IPick and Vilmerl 
119871 [20081 120081 for recent reviews). Basic physical understanding of the gener- 
ation of type III radio burst EM waves in plasma by an electron beam has been 
with us for over five decades (Ginzburg and Zhclczniakov, 1958) and involves 
the generation of Langmuir waves by bump-on-tail unstable electron distribu- 
tions and subsequent mode conversion of the longitudinal Langmuir waves into 
escaping, transverse EM radiation. Theoretical efforts in understanding of type 
III solar radio bursts can be grouped into three categories: 

(1) Quasilinear theory of type III sources that use kinetic Fokker-Planck type 
equation for describing the dynamics of an electron beam, coupled with spectral 
energy density evolutionary equations for Langmuir and ion-sound waves have 
long been studied. The most essential result is that the spectral energy density 
of the Langmuir wave packets (that are excited by the bump-on-tail unstable 
beam) travels along the open magnetic field lines with a constant speed and this 
is despite the quasilinear relaxation (formation of a plateau in the longitudinal, 
along the beam injection direction, velocity phase-space of the electron distri- 
bution function), hence, this implies some sort of beam marginal stabilisation 
QKaplan and Tsytovich, 1968||Smith, 1970[|Zaitsev, Mityakov, and Rapoport, 1972| | 



Mel'Nik, Lapshin, and Kontar, 1999[|Kontar and Pecsc li, 2002). Inclusion of EM| 



emission component into the quasilinear theory in some models is based on so- 
called drift approximation ( [Hillaris, Alissandrakis, and Vlahos, 1988 ; Hillari s~e£ al., 1990} [ 
|Hillaris et al., 1999[ ), where nonlinear beam stabilisation during its propagation 
(so called free streaming) is based on Langmuir-ion acoustic wave coupling via 
pondcr-motivc force and EM emission is prescribed by a power law of the 
beam to ambient plasma number density ratio. These models can be success- 
fully compared with the observed dynamical spectra to constrain key model 
parameters. 

(2) Stochastic growth theory ( |Robinson, 1992 Robinson, Cairns, and Gurnett, 1992| ),[ 



in which density irregularities induce random growth, such that Langmuir waves 
are generated stochastically and quasilinear interactions within these Langmuir 
clumps cause the beam to fluctuate about marginal stability. Further this ap- 
proach has been developed into a numerical simulation tool that effectively can 
reproduce observational features of the type III bursts ( |Li, Cairns, and Robinson, 2008[ ) .| 

(3) Full kinetic simulation approach of type III bursts ( |Kasaba, Matsumoto, and Omura, 2001| [ 



|Sakai, Kitamoto, and Saito, 2005| |Rhee et al., 2009| |Umeda, 20101 ) to this date 
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used Particle-In-Cell (PIC) numerical method. This effort is mainly focused on 
understanding of basic physics rather than direct comparison with the observa- 
tions because the size of simulation domain of the models corresponds to only 
few 1000 Debye lengths which is roughly l/10 10 th of 1 AU. Thus, such models 
deal with micro-scales. 

Each of the above theoretical approaches have their advantages and disad- 
vantages. For example as it is well explained in ICairnsl (|1985p . in quasilinear 
theory the particle distribution function is split into a slowly varying average 
part and a rapidly varying part due to a wave. Then major approximation is 
that interactions between wave modes are neglected and only back-reaction of 
the waves on the slowly varying average part of the distribution function is con- 
sidered. Therefore, if computational resources permit, the full kinetic approach 
is more desirable. However, to this date only PIC method was used. It is well 
known that PIC approach suffers from large, so called, "shot-noise" level of which 
scales as one over the root of number of particles. Moreover, in PIC approach 
typically there are few hundred particles per cell (which is normally one Debye 
length long in fully electromagnetic PIC codes). In Vlasov-Maxwell approach 
instead of solving for individual particle dynamics, without loss of generality 
(or any kinetic physics), collisionless Vlasov equation for electrons and ions is 
solved in which EM fields are self-consistent. In this study we use 80x80 velocity 
grid which in PIC equivalent would be having 80 x 80 = 6400 particles per cell 
(instead of few hundred) . We have also done convergence tests by increasing the 
resolution both in velocity space and spatial domain and confirmed the results' 
convergence. This shows that Vlasov-Maxwell approach probes more finely phase 
space of the problem but this comes at a substantial memory cost. 

The paper is organised as following. In Section 2 we present the model and 
main results, including (i) a numerical run of the inhomogeneous plasma without 
a beam which turns out to be aperiodically unstable to a Larmor drift. The latter 
results in Langmuir and EM wave generation, as well as density filamcntation 
(Section 2.1); (ii) a numerical run of the homogeneous plasma with a low density 
beam which generates Langmuir and EM wave via plasma emission mechanism 
(Section 2.2); and a numerical run of the inhomogeneous plasma with a high 
density beam in order to combine the effects of both density inhomogencity and 
the presence of the beam (Section 2.3). 

2. The model and general theoretical considerations 

Our numerical model is implemented using a rclativistic, fully electromagnetic 
Vlasov-Maxwell code called VALIS ( |Sircombe and Arber, 2009[ ). The code is 
using a conservative, split-Eulcrian scheme based on the Piecewise Parabolic 
Method for the update of the particle distribution function and utilising the exact 
particle fluxes to calculate the current in the solution of Maxwell's equations. 
In particular relativistic Vlasov's equation is solved for species a (a = e,i for 
electrons and ions respectively) 




(1) 
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in conjunction with the Maxwell's equations 



V-E = (0 / £o , V-B = 0, VxE = - — , VxB = ,i„J+ T - (2) 



where charge and current densities are specified in self-consistent manner, 



VALIS is 2D2V code in that it has two spatial dimensions {x, y) and two 
corresponding velocity components (u x ,u y ), while electric and magnetic field 
components that are solved for are (E x ,E y , 0) and (0, 0, B z ). Distance and time 
are normalised to c/ui pe and u>~ e , while electric and magnetic fields to uj pe cm e /e 
and oj pe m e /e respectively. Temperature is normalised to to m e c 2 /k. Here ui pe = 

n e e 2 / (eom e ) is the electron plasma frequency, n a = J f a d 3 u is the number 
density and all other symbols have their usual meaning. 

We intend to consider a single plasma thread (i.e. to use 1.5D geometry), 
therefore space components considered are (x,y) — (25000Ad, IAd) with A^> = 
v t h,e/u pe being Debye length (here Vth.e = y/kT/m e is electron thermal speed). 
We would like to resolve full plasma kinetics, therefore we set spatial grid size 
as IXd- In practice this means we set plasma temperature at T = 10 5 K (i.e. 
fix v t h.e at 4.12 x 10 _3 c), which corresponds to high solar corona, above active 
regions. In the presented results (and in the VALIS code generally) spatial scales 
are normalised to c/u pe . We do not fix plasma number density and hence uj pe 
deliberately, because we wish our results to stay general. In order to achieve this 
generality (and consistency of the results) it is important to keep normalised 
B z o = B z q/ {uj pe m e / e) = 0.01 the same in all numerical runs. Because B z q is 
normalised to w pe m e /e, no matter how plasma density and hence uj pe changes, 
(i) ratio of Debye length and c/uj pe , i.e. Xd/ (c/uj pe ) = 4.12 x 10~ 3 and ratio 
of electron Larmor radius and c/ui pe , i.e. r^^/ (c/ui pe ) = 4.12 x 10 _1 stay the 
same. Such choice means, of course, that magnetic field in Tesla is variable. For 
example, if we set plasma number density to hq — 10 15 m™ 3 (i.e. fix ui pe = 1.78 x 
10 9 Hz radian), this sets Debye length at A^ = 6.90 x 10~ 4 m=4.12 x 10" 3 c/wp e 
and electron Larmor radius at r/,. e = 6.90 x 10~ 2 m=4.12 x 10 _1 c/w pe . Also, 
then B z q = 1.01 x 10~ 4 T ps 1 gauss. If we set plasma number density to no = 
10 -5 m -3 5 this scts De b ye length at X D = 6.90 x 10 6 m=4.12 x 10- 3 c/w pe and 
electron Larmor radius at rx je = 6.90 x 10 8 m=4.12 x 10 _1 c/w pe - Also, then 
B z q = 1.01 x 10 -14 Tesla. In other words, appropriately adjusting plasma number 
density no, physical domain can have arbitrary size e.g. Sun-earth distance (but 
then unrealistically low density has to be assumed). Since the number of grid 
points and domain size (normalised to c/ui pe ) is set independently, we have to 
make sure that grid size is IAd by setting n x = 25000 and L x ^ nax = 25000x A_d = 
102.94c/w pe and n y = 1 and L VtTaax = 1 x Xn = 4.12 x 10~ 3 c/o;p e . For the 
velocity we have (u x , u v ) = (80, 80) grid points with maximal possible velocities 
for electrons allowed set u x , max = u y , max = 0.25c (u x , m ax = u y ^ max = 0.4c in 
Section 2.3) for electrons and u x . max = u y ^ max = 0.25/vT836c = 5.83 x 10 _3 c 




(3) 



and 7* = ^1 + |u| 2 /c 2 . 
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for ions. Since the code does not allow to set magnetic field along x-axis (because 
B = (0,0, B z )), to represent the situation adequately, we have only a choice to 
set a B z q component, which we fix at 0.01 in normalised units. This is not an 
unreasonable value for, transverse to the considered plasma thread, component 
of magnetic field above an active region. At first, it seems unrealistic to ignore 
magnetic field along x. However, firstly, bulk of the work in the quasilinear 
theory indeed makes the same assumption (ignores longitudinal magnetic field). 
Secondly, it is known ( |Alexandrov, Bogclan kevich, and Rukhadze, 1988[ ) that in 
the case of weak fields, general picture of excitation of the Langmuir waves 
by a low density (rib ^ n e) electron beam (i.e. their resonant interaction) via 
Chcrenkov resonance is not much different from the case without the magnetic 
field. 

In all presented numerical runs boundary conditions for all quantities are 
periodic. A typical numerical run takes 32 hours on 256 processor cores (Dual 
Quad-core Xeon, eight cores per computing node) . 

We now briefly re-iterate key facts about beam-plasma interaction theoretical 
framework. In the case without magnetic field (or in the weak field case) cold 

plasma dispersion relation yields two possible modes (e.g. |Alexandrov , Bogda nkevich, and Rukhadze| [ 
115551 p.156): 

tJ 2 = k 2 c 2 +u 2 +oj 2 pb 7- 1 , (4) 



k 2 ± v 2 j 21 



UJ 2 (u> — fc||Wf,) 2 



(5) 



Here, cu p b = \Jribe 2 j (eom e ) the beam plasma frequency and Vb is its speed. 7 is 
the usual Lorcntz factor for the bulk motion of the beam (in cold plasma ap- 
proximation thermal motions of plasma are absent). Note that here only electron 
and beam contributions are retained whilst ion contribution is ignored due to its 
smallness. Equation (4) describes a stable, purely transverse (Elk), EM wave 
which does not interact with the beam because E • vj, = 0. If Bo \\ z then in 
the cold plasma approximation beam can only propagate along z-axis (Note that 
our numerical simulations are with finite temperature). EM wave described by 
Equation (4) has non zero E y component only. Equation (5) describes oblique 
wave which has both E^ = E z and E x components and hence can interact with 
the beam via Cherenkov resonance. By putting uj ~ k^Vb+S = k z Vb+S into Equa- 
tion (5) growth rates (for k^Vb < uj pe when the oblique mode becomes unstable), 
<5, can be easily found: (i) away from the plasma frequency, u> 2 w k 2 v 2 ^ uj 2 e 
(non- resonant case), 



-3/2 



"fcjf + k\ 1 2 



1/2 



^l-u 2 J(k 2 v 2 ) 



Im(6 nr ) CX LU pe 



lib 



1/2 



(6) 
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(ii) close to the plasma frequency, ur ps kiv£ w Wp e (resonant case), 

1/3 



5 r 



1 ^ + fci7 2n1 ' 



2n e 7 3 A: 



Im(S r ) oc w pe ( ^ ) . (7) 



Naturally, the resonant growth rate is much larger than the non-resonant one, 
S r ^> S nr , due to the beam's low density (n& <C n e ). As is also clear from Equation 
(5), the Cherenkov resonance effectively excites Langmuir (longitudinal) waves 
with the dispersion relation cu 2 ps cu 2 e (note that thermal effects are ignored in 
Equations (4)- (9)). 

In the case of strong magnetic field (u ce 3> w pe ) there are two types of waves, 
dispersion relations of which are given by ([Alcxandrov, Bogdank cvich, and Rukhadze, 1988] )[ 

k 2 - oj 2 /c 2 = 0, (8) 

k\ + [k\ u 2 /c 2 ){l u 2 Ju 2 - u 2 pbl -i/{u k\\v b ) 2 ) = 0. (9) 

Our Equation (9) without the beam contribution term is also derived bv lArons and Barnard! ] 

119861 their Equation (49), and is referred to as O-mode. The beam does not 

interact with the purely transverse EM wave given by Equation (8), while it can 

interact with the slow wave w_ given by Equation (9) (Equation (9) describes 

oblique EM waves which have both E\\ and E x and in the absence of the beam 

reduces to the fast (w+) and slow (w_) modes with dispersion 

Electron beam then can 



2 pe + k 2 c 2 ± J{u 2 pe + k 2 c 2 ) 2 - ico^kl 



c- 



resonantly interact with the slow mode when w_ intersects with u> = kuVb line. 
Thus, the presence of the longitudinal magnetic field can only alter value of kit at 
which the Cherenkov resonance occurs. It should be stressed that Equations (4)- 
(9) are derived in the case electron beam propagating strictly along the magnetic 
field. Therefore, beam is strictly decoupled from the purely transverse EM wave 
(Equation (8)). If the beam has a small k± (note that k± in Equations (4-9) refers 
to that of a wave mode) then it can couple to EM wave and hence generate it. 
This is our main motivation for have small B z q so that when beam is injected 
along x-axis it can couple to the EM mode. Also, Hsu (2010) showed, using 
relativistic Vlasov equation, that in unmagnetised plasma, EM and plasma wave 
conversion efficiency diminishes to zero at both 0° and 90° incidence angles and 
peaks between 10 — 20° depending on plasma temperature. Here degrees refer to 
angle between wavenumber k and density gradient direction. Resuming aforesaid 
we acknowledge that the present numerical model is not suitable for describing 
the type III solar bursts directly. Ideally, to represent the true physical reality, 
it would be preferable to set large B x q in addition to small B z q (as in solar wind 
Parker spiral in the upper solar corona). However, since the VALIS code only 
solves for (E x ,E y ,0) and (0, 0, B z ) we have to simply make sure that when the 
beam is injected along cc-axis it can couple to the EM mode in order to capture 
the essential physics. At the same time we stress that existence of B x q is not a 
requirement for the generation of type III bursts per se, what is essential is to 
have finite k± in the beam so that it couples to the EM wave (here we achieve 
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this by setting small B z q only) . We also note that it was our intension to consider 
magnetised plasma with the beam injected strictly along the physical domain. In 
principle, the coupling to EM wave could have been also achieved by switching 
off the magnetic field altogether and in addition to uq x , setting uq v = 0.2 — 0.5c. 
This would have created a situation with non-zero k± too, thus facilitating the 
coupling of the beam to an EM wave. 

VALIS code allows to set any desired number of plasma particle species. 
Therefore, because we intend to study spatially localised electron beam on top of 
the inhomogeneous or homogeneous Maxwellian electron-ion plasma, we solve for 
three plasma species electrons, ions and the electron beam. The dynamics of the 
three species, which all interact via EM interaction, can be tracked independently 
in the numerical code. Velocity distribution function for electrons and ions is 
always set to 

f e Au x ,u y ) = e- m ^<+<^ 2T \ (10) 

where m rie — 1 for electrons and m r ,i = 1836 for ions. When cases with the 
beam are considered we set the following distribution 

f b (u x ,u y ) = h b e-^-°^ 2+ <y^. (11) 

where fib is normalised beam number density (fib = nb/n e Q ) and it is fib = 
5 x 10~ 6 for low density beam (Section 2.2) and fib = 5 x 10~ 2 for the dense 
beam (Section 2.3). The normalised number density of the background plasma 
in the homogeneous case (Section 2.2) is no = 1. Thermal spread of the electron 
beam is specified by setting Tb = 9T = 9.0 x 10 5 K. Note that the beam is injected 
along the x-axis, transverse to the background magnetic field B z q. Physics of the 
initiation of the beam is believed to be related to the magnetic reconnection. In 
2D case reconnection electric field at a magnetic null is in direction out-of-plane 
where magnetic field lies. Therefore, it is not unreasonable to consider situation 
when beam is injected as in our model. Also, beam injection transverse to the 
magnetic field can result from accelerated electrons from the collapsing magnetic 
traps (jKarlicky and Kosugi, 2004). In the inhomogeneous cases (Sections 2.1 and 



2.3) background plasma normalised number density is set to 



n (x) = 1/ 



1 + 10 8 e 



-[(*-L x , max /2)/21] 4 



(12) 



This density profile mimics a factor of 10 8 density drop from the corona no = 
10 15 m -3 to tiau = 10 7 m -3 at 1 AU. Because it is known that numerically 
most precisely implemcntablc boundary conditions are periodic ones this density 
profile effectively represents mirror-periodic situation when the domain size is 
doubled, i.e. at n (x = 0) = n (x = L x , max ) = 1 while n (x = L x , max /2) = 
10 -8 . This way "useful" or "working" part of the simulation domain is < 
x < L x ^ max /2. Spatial width of the density gradient is Lm ~ 5c/w pe (see e.g. 
Figure 7(c), thick solid curve for 90 < x < 95). When cases with the beam are 
considered we set its following density profile: 

n b (x) = n 6 e-^- 5 '/ 3 ] 4 (13) 
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which means that the beam is injected at x = 5c/u!p e and its full width at half 
maximum (FWHM) is also ps 5c/ui pe (see Figure 5(c) dotted curve). Plasma 
beta in this study, based on the above parameters, is set to fi = c 2 /v\ = 
V thi/ V A = ( l 'th,i/c) 2 (u! p i/uJci) 2 = 0.17. (c s and va are sound and Alfven speeds 
respectively.) It should be noted that, strictly speaking, pressure balance in the 
initial conditions is not kept. There are two reasons for this: (i) solar wind is 
not in "pressure balance" and it is a continually expanding solar atmosphere 
solution; (ii) plasma beta is small therefore it is not crucial to keep thermody- 
namic pressure in balance (because its effect on total balance is negligible) and 
the initial background density stays intact throughout the simulation time (see 
e.g. Figure 7(c), thick solid curve for 90 < x < 95). 

2.1. Larmor drift- unstable case, inhomogencous plasma without a beam 

It is well known that the mode conversion from electrostatic to EM waves 
near the plasma frequency is possible by linear coupling on a density gradient. 
lYin et a~L\ (|1998[) examined the mode conversion from electrostatic to EM waves 
near the plasma frequency in the Earth's electron fore-shock. The conversion 
and reflection coefficients were obtained by solving coupled differential equa- 
tions in a weakly magnetised warm plasma with a longitudinal linear density 
gradient. Results indicated that the fore-shock first harmonic EM emissions 
and the backward-propagating Langmuir waves required for the generation of 
the second harmonic EM waves could be efficiently generated by the linear 
conversion process in an inhomogeneous plasma. Therefore, originally the aim 
was to study super-thermal beam injection into plasma with homogeneous and 
inhomogeneous plasma to study the effect of the density gradient on the level of 
EM wave generation. However, we found originally unforeseen outcome in that 
with or without electron beam background density gradient regions generate 
perturbations in all quantities n(x), E x , E y and B z . The results are presented 
in time-distance plots in Figure (1) and Figures (2)-(3). Dynamical picture is 
presented in movie 1 in the electronic supplement to this article. 

As shown below, the obtained results can be interpreted by means of Larmor- 
drift unstable mode Alexandrov, Bogdankevich, and Rukhadze[[1988l p. 239. There-| 



fore before discussing this numerical run results, let us briefly summarise key 
facts about the Larmor-drift instability. For frequencies smaller that the Larmor- 
drift frequency (J < ujld = k y v 2 h a /{ui ca LiH) this mode is aperiodically unstable 
in certain regimes (which as we will see below are always taking place in our 
model). Physical meaning of Larmor-drift in inhomogencous plasma is clear. 
When magnetic field is directed along z and plasma inhomogeneity is along x- 
axis, variation of the particle Larmor radii (due to the inhomogeneity) generates 
transverse to the both directions current J y ps q a n a v 2 h a j 1 (lu COi Lih)- This drift 
is not related to the actual motion of centres of the Larmor orbits, and it is 
preferentially realised in low beta plasmas when particles are magnetised. It is 
important to note that such Larmor-drift instabilities may occur in Maxwellian 
plasmas and they are not related to the existence of a positive-sloped region in 
the velocity distribution function. In this sense the instability can be regarded 
as hydrodynamic. In the limit of long wavelength approximation, Ax rt,%i 
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Figure 3. (a) E x (x,t = 50), (b) E y (x,t = 40), (c) n e {x,t = 50)-n e0 , (d) _B z (x,i = 40)-B z0 , 
(e) rii(x,t = 50) — n^o, (f) fe(v x )- In (f) dotted curve represents f e (v x ,t = 0) while solid one 
fe(v x ,t = 50). For E v and £? z — B z o the time snapshot t = 40 is chosen before the end of the 
simulation time, t = 50, in order to show the spatial profiles before EM fronts collide. It can 
be seen in Figure 1(c) (look horizontally across t 50) that EM fronts collide at about t Rj 50. 



and oj 3> k z v t h,e, (which naturally holds because in our case k z — > because 
our domain size is infinite in z-dircction); also when u> <C c^ld', w and 
ujpi 3> ui c i (for our choice of parameters uj p i/ui C i — c/va = 4.29 x 10 3 >> 1 and 
this holds for arbitrary background plasma number density) . IRukhadze and Silinl 
(|1964p derived following dispersion relation for the Larmor-drift mode 

„ 2 = „2 k z T e m P dlnn e T e 
c% k\ Ti m e dlnn.iTi 

Here dkiA/dlnB = (d\nA/dx)/{d\nB/dx) = {BdA/dx)/(AdB /dx) notation 
is used. It is clear that uj <C ^ld and to <C uj C i conditions hold too (confirming in 
retrospect) because k z — > 0. Equation (14) shows that condition for instability 
is 

§S§>° < 15 > 

which for uniform temperature plasma T e = Tj = const reduces to 
d In n e /d In = 1 > is always fulfilled when density variation for ions and 
electrons is the same. Note that the plasma will be drift unstable for both 
increasing (positive) and decreasing (negative) density gradients because the ra- 
tio d\\\n e / dlnrii will be always positive (negative/negative or positive/positive 
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is positive). So, we conclude that our inhomogeneous plasma set up is always 

Larmor-drift, aperiodically unstable. It is also interesting to note that ( |Alexandrov, Bo gdankevich, and Rukha 
p. 169 showed that aperiodic instabilities can lead to density filamentation (cre- 
ating of spatially thin threads) . We can indeed see similar filamentary structures 
in density (and E x ) in Figures 1(a), 1(b) and 3(a), 3(c). 

We gather from Figures 1(a) and 1(b), that E x and n e — n e Q perturbations 
travel rather slowly compared to E y and the fast part of B z — B z q (dark oblique 
strips in Figure 1(d)). As can be inferred from Equation (12) and thick solid 
curve in Figure 7(c) (for 90 < x < 95), the lengthscale of the background density 
gradient is Ljh ~ 5c/Lu pe which roughly corresponds to the distance travelled 
by E x and n e — n e0 perturbations (see e.g. start and end positions of rightmost 
bright streak in Figure 1(a) or location of the rightmost peaks in Figures 3(a) 
and 3(c)), i.e. 11 — 6 = 5c/uj pe . This distance is travelled in time 50uj pi }. Thus, 
the phase speed is estimated as 0.1c. Generally E x and n e — n e Q perturbations 
repeat the shape of the density gradient, and other runs (not shown here), with 
varied density gradient strength, confirm this. E y perturbation as can be seen 
from Figure 1(c) travels from the density gradient edges with a speed « c (as 
the slope of dark and bright bands is unity). B z — B z0 (Figures 1(d) and 3(d)) 
perturbation has two parts: slow moving part that travels with speed 0.1c (as in 
E x and n e — n e o) and and smaller, leading pulses which travel with speed of c. 
We also gather from 3(f) that electron temperature is slightly increased (i.e. the 
electron distribution function gets broader at t = 50^^ (solid curve) compared 
to t = (dotted curve)). 

To estimate frequency both of the slow and fast perturbations, we note the 
number of bright features along e.g. left edge at x = 5 in Figure 1(a) which is 7, 
counting from the first. Note that the density gradient left edge where no w 1 is 
at x = 5. The frequency estimate can be better done using Figure 2 where we 
plot E x [x = 5, i) and E y (x = 5, t). The estimate is as follows. In Figure 2(a) the 
time difference between leftmost and rightmost peaks is At = 47.2—3 = 44.2w pe 1 . 
Thus 7(1//) = 44.2W- 1 = 44.2 x (2irf pe )- 1 , and / ~ f pe . Similar calculation 
for Figure 2(b) yields the time difference between leftmost and rightmost peaks 
is At = 47.5 - 5 = 42.5a;- 1 . Thus 7(1//) = 42. 5W" 1 = 42.5 x (2 7 r/j, e )" 1 , 
and / = 1.035/ pe . We therefore conclude that E x oscillates at local plasma 
frequency and corresponds to a plasma wave. Whereas E y perturbations are EM 
type (escaping radiation) and oscillate just above the plasma frequency 1.035/ pe . 

We would like to stress that the generation of perturbations in the considered 
Larmor drift unstable case is not due to the fact that pressure balance is not kept. 
In our case plasma beta is small, therefore it is not crucial to keep thermodynamic 
pressure in balance and the initial background density stays intact throughout 
the simulation time (see e.g. Figure 7(c), thick solid curve for 90 < x < 95). We 
have performed numerical runs where temperature was varied as inverse of no (x) 
so that po = no(x)kBTo(x) = const (B z q is constant throughout this study) 
yielding perfect total pressure balance - similar approach to pressure balance 
was adopted by |Tsiklauri, Sakai, and Saito (|2005|) . We confirm that even when 



total pressure balance was kept, Larmor drift instability still developed and 
physical system behaviour was similar to what is presented here. 
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Figure 4. 
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As in Figure 1 but for the case of homogeneous plasma with low density beam. 



We note that if finite Bq v = Bq± is added, this will change the Larmor-drift 
stability criterion. This instability will stabilise if (Rukhadzc a nd Silin, 1964[ ): 



B 



> 



Vth, 



LOriL 



ci^IH 



(16) 



We plan to study this stabilisation issue in a following publication (work in 
progress), using EPOCH 1.5D particle-in-cell code which allows to choose all 
background magnetic field components (not as VALIS, which only allows to 
consider (E x ,E y ,0) and (0,0, B z )). 

We conclude this subsection with the observation that we found a new pos- 
sibility for exciting plasma frequency EM radiation by means of a universal, 
aperiodic Larmor-drift instability. By universal we mean that satisfying the 
instability criterion (see Equation (14)) is quite plausible in many astrophysical 
situations. Condition lu <C w C i requires that k 2 7 Jk\ ~ L\jL 2 z -C m e /m p i.e. 
length of the domain should be at least ~ 43 times longer than its width. 

2.2. Plasma emission case, homogeneous plasma with low density beam 

To suppress the Larmor-drift instability (because we cannot achieve this by 
imposing suitable Bq x ) we now set uniform normalised plasma number density 
to no = 1 and inject a low density beam with the parameters specified above. 

The results of this numerical run are presented in Figures 4 and 5, while time 
dynamics is presented in movie 2 (see the accompanying electronic supplemen- 



SOLA: cesra2010_dt.tex; 21 October 2010; 0:45; p. 12 



Vlasov-Maxwell simulations of radio emission in solar corona 




-0.2 -0.1 0.0 0.1 0.2 -0.002 -0.001 0.000 0.001 0.002 



Figure 5. (a) E x (x,t = 50), (b) E y (x,t = 50), (c) n e (x,t = 50) (thick solid horizontal line), 
ni,(x, t = 0) (dotted curve) and nt,(x, t = 50) (thin solid curve), (note that n;, was scaled by a 
factor of 2 X 10 5 to make it visible), (d) B z (x,t = 50) — B z o, ( e ) fe(v x ,t = 0) (dotted curve) 
and f e (v x ,t = 50) (solid curve), (f) fi(v x ,t = 50) (solid curve) (fi(v x ,t = 0) is also plotted 
with a dotted curve, but to a plotting accuracy the curves overlap, indicating no ion heating 
takes place). 

tary material). One striking novel feature immediately seen in Figures 4(a) and 
4(b) (E x and n e — n e Q where n's include initially injected beam contribution) is 
the excitation of standing waves in the spatial location of the beam injection. By 
counting the number of bright features in the elapsed time, it is clear that the 
oscillations are at about Lu pe . Moreover, E x oscillates as one solid feature in the 
spatial location of beam injection (oscillation spatial width coincides with the 
beam width). While n e — n e Q also exhibits standing waves, but these are in the 
regions of positive and negative density gradients of the back and front of the 
beam (recall that here background plasma is homogeneous). These oscillations 
are in anti-phase, i.e. at given t = const over-density and under-density is ob- 
served. As with Larmor-drift instability (Section 2.1) this was unforeseen result. 
However, again literature analysis enabled us to find a suitable interpretation. 
There are two possible regimes. If the phase of the waves is locked, in the strong 
instability regime, waves can appear with the frequency close to uj pe in the 
location of the beam injection (Pavlcn ko and Petviashvili, 1977[ ). In the case of 
a weak turbulence regime, the formation of strong Langmuir waves also observed 
near the beam injection site ( |MerNik, Lapshin, and Kontar, 1999[ ) (according to 
their Equation (15) Langmuir turbulence spectral energy density near the beam 
injection point depends on the phase velocity as W oc v 5 ). In the electromagnetic 
E y and B z — B z q components (Figures 4(c) and 4(d)) again standing wave at 
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the beam injection location can be also seen, but in addition this serves as a 
source to the escaping EM radiation. These can be seen as oblique dark and 
bright strips with a slope close to unity, thus propagating at about speed of 
light. Wide dark oblique strip (with narrow bright front) in Figure 4(d) with a 
slope Ax/At = (15 — 5)/50 = 0.2 corresponds to the beam wake (recall that 
beam travels with speed 0.2c). Note that in 4(c) and 4(d) there are also EM 
waves present near x ~ 100. This is because the standing wave centred on 
x = 5 generates EM waves travelling in both directions. Because of the periodic 
boundary conditions, waves that travel to the left, appear on the right side of 
the simulation domain. 

Figure 5 provides further details: 5(a) shows E x at time t = and is 

made of two parts (i) a deep centred on x = 5 corresponds to the standing 
wave at plasma frequency and it nicely coincides with beam injection site, see 
dotted pulse in 5(c) which shows the beam at t = 0; (ii) a smaller hump in 
5(a) at x w 16 which roughly coincides with the location to where beam has 
travelled in time t = 500;"/ (solid pulse in Figure 5(c) centred at x rj 15). 
Figure 5(e) indicates that bulk plasma electron distribution function heats up 
(solid peak centred on v e — which corresponds to t = bQui~^ is wider than 
at t = 0). Also we can see that the beam slows down from 0.2c to 0.17c (small 
bump (dotted curve) centred at 0.2c which represents the beam at t = shifts 
to 0.17c (solid curve bump) which is the beam at t = 50w pe 1 ). We note that 
since the beam is nine times hotter than the background plasma, Tj = 9T = 
9.0 x 10 5 K, Larmor radius of the beam is three times larger than that of back- 
ground plasma, ri,,b = 1.24c/w pe . This is smaller than the distance traveled by 
the beam (ss 10c/w pe cf. Figure 5(c)). Thus, the beam partial magnetisation 
can be regarded as the main cause of its slowing down. We also observe that 
there is no substantial quasilinear relaxation (i.e. plateau does not form) which 
corroborates basic features of the quasilinear theory. This is due to the fact that 
the growth rate of resonant Langmuir waves given by Equation (7) is small, 
as in this run nf,/?i e = 5 x 10~ 6 . A simpler estimate for quasilinear relaxation 
time, r, (time of establishing the plateau) is achieved by using r = n e / (nbU> pe ) 
(e.g. ( jMel'Nik, Lapshin, and Kontar, 1999[ )). Based on this, we see that in our 
case r = 2 x lO 5 ^^ 1 . Thus, it is not surprising that in 50 plasma frequencies 



we do not see substantial quasilinear relaxation. Mel'Nik, Lapshin, and Kontar 
(1999]) quote the criterion of weak turbulence regime of quasilinear theory to 
apply as e = ntm e vl / (nom e vf h e ) <C 1. Here, e « 10~ 2 <C 1, thus we are well 
in the quasilinear regime. Another interesting corroboration of the quasilinear- 
thcory is that shape of the beam does not change i.e. despite the fact that 
small density beam plunges through the background plasma at a speed of 0.2c, 
it stays intact. |Mel'Nik, Lapshin, and Kontar| (|1999[) offer suitable explanation 
for this fact based on the beam particle kinematics. To avoid duplication we 
refer the interested reader to this paper for the details. Figure 5(f) confirms that 
despite the fact ions are treated in the numerical code as moving, still there is 
no significant change in their velocity distribution function. 

It is interesting to note that newly established standing waves can offer an 
alternative explanation for the horizontal strips observed in some dynamical 
spectra. lAurass ~ei al. (2010) report a narrow-band, short duration line emission 
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at 314 MHz which is interpreted as a gyro-resonance line. We note that the 
observed feature can also be explained by an EM emission emanating from 
the standing waves. This naturally explains the fact that there is no drift in 
frequency, as the standing wave remains in the same spatial location, hence 
there is no change in density and in emission (plasma) frequency Moreover, 
if we look at the dynamical spectra from Figure 1 from lAurass et aU PUP]) 
we see that initially, the line intensity increases in time. This behaviour is 
very similar to what is seen our 4(c) where the intensity increase at x — 5 
in time can be seen. We can roughly estimate intensity of the line predicted 
by our model as follows. From Figures 5(b) and 5(d) we can read off typi- 
cal amplitudes of the escaping EM radiation in the form of standing wave as 
E y ~ 1CP 7 and B' z = B z — B z0 ~ 2 x 10~ 7 . Flux of the Poynting vector is than 
F = |F| ~ E y B z /[j,Q. Recovering physical units from the normalised quantities 
yields F«5x 1CP 4 Wm~ 2 . Assuming the area of the emitting source is A =20 
Mmxl Mm and plasma frequency of f p — 300 MHz, we can estimate the flux 
density at 1 AU as F d = FA/[^(lAU) 2 f p ] « 1.2 x 10" 22 Wm~ 2 Hz" 1 (w 1 
sfu) . This is of the order of the lAurass et al.\ (|2010[) estimate of few sfu for the 
line flux density. 

In summary, the plasma emission case in a homogeneous plasma with low 
density beam confirms the general picture of quasilinear theory even in the case 
when the beam is injected transversely to the magnetic field, with a main novelty 
being established that the spatial location where beam is injected serves as a 
source for standing plasma waves which, in turn, generate escaping EM radiation 
that oscillats at plasma frequency. This offers a possible new interpretation of 
the horizontal strips observed in some dynamical spectra. 

2.3. Larmor drift-unstable, plasma emission case, inhomogencous plasma with 
dense beam 

We now combine both effects considered in Sections 2.1 and 2.2 by considering 
Larmor drift-unstable plasma emission case in which plasma is inhomogeneous 
and dense beam is injected transverse to the magnetic field and along the density 
gradient. 

We gather from Figure 6(a) that E x (plasma wave component) now comprises 
of two parts: a weak standing wave at the location of the beam injection which 
oscillates with frequency uj pe and a strong wake created by the beam which now 
shows much more dispersion, as we depart from the quasilinear regime where 
non-linear interactions between wave modes are ignored. Time-distance plot for 
n e — n e0 (Figure 6(b) is dominated by a wake created by the beam). EM wave 
components (E y and B z — B z q) as in Section 2.2 show similar behaviour where 
the standing wave centred on x = 5 generates escaping EM waves travelling in 
both directions. 

Figure 7(a) shows that E x has two parts: a small leftmost bump which is 
from a standing plasma wave and a beam generated perturbation which travels 
roughly with the beam speed. Figures 7(b) and 7(d) for E y and B z — B zQ , show 
in more detail, how the standing wave centred on x = 5 generates EM waves 
travelling in both directions with speed of light. We gather from Figure 7(c) that 
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as the dense beam plunges into the plasma, it no longer retains its shape as in 
quasilinear regime (compare to Figure 5(c)). Figures 7(e) and 7(f) show dynamics 
of electron and ion distribution functions. We see that as the beam/background 
plasma number density ratio, rib/n e = 5 x 10~ 2 , is no longer small, two effects 
can be observed: (i) quasilinear relaxation happens very fast; and (ii) substantial 
electron return current (wide wing with negative velocities in Figure 7(e)) is 
generated. For ions sizable heating takes place (Figure 7(f)). Dynamical picture 
of the system evolution can be studied using movie 3 (see the accompanying 
electronic supplementary material). In the considered case r = 20w~ e 1 . Thus, it 
is not surprising that in 50 lo~} we see substantial quasilinear relaxation taking 
place. The criterion of weak turbulence regime e = ni,m e v^ / (riom e v^ h ) <C 1 is 
no longer fulfilled as here e « 10 2 3> 1. Thus the physical system is not in the 
quasilinear regime. 

In summary, we see for this set of results that the system is driven by the 
effects of the beam while Larmor-drift effect, whilst present, is not dominant. 
Significant deviation from the quasilinear theory is found which manifests itself in 
(i) fast quasilinear relaxation, (ii) disintegration of the beam, and (iii) generation 
of significant electron return current and ion heating. 
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Figure 7. (a) E x (x,t = 50), (b) E y (x,t = 50), (c) n e (x,t = 50) (thick solid curve), 
ni,(x,t = 0) (dotted curve) and ni,(x,t = 50) (thin solid curve), (note that n;, was scaled 
by a factor of 20 to make it visible), (d) B z (x,t = 50) — B z q, (e) f £ (v x ,t = 0) (dotted curve) 
and f e (v x ,t = 50) (solid curve), (f) fi(v x ,t = 0) (dotted curve) and fi(v x ,t = 50) (solid 
curve) . 



3. Conclusions 

We used 1.5D Vlasov-Maxwell simulations to model EM emission generation 
in a fully kinetic model for the first time in the solar physics context. We 
studied plasma emission mechanism and Larmor drift instability in a single 
plasma thread that joins the Sun to Earth with the spatial scales compressed 
appropriately. The results can be summarised in three points: 

• We established that 1.5D inhomogeneous plasma with a uniform back- 
ground magnetic field directed transverse to the density gradient is ape- 
riodically unstable to Larmor-drift instability. This instability results in 
a novel effect of generation of EM emission at plasma frequency. The 
generated perturbations consist of two parts: (i) non-escaping (trapped) 
Langmuir type oscillations which are are localised in the regions of density 
inhomogencity, are highly filamentary, and the period of appearance of the 
filaments is close to electron plasma frequency in the dense regions; and (ii) 
escaping electromagnetic radiation with phase speeds close to the speed of 
light. 

• In the uniform density plasma case (when plasma becomes stable to Larmor 
drift instability), when a low density super-thermal, hot beam is injected 
along the domain, in the direction perpendicular to the magnetic field, 
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plasma emission mechanism generates non-escaping Langmuir type stand- 
ing oscillations at the beam injection location, which in turn generate 
escaping electromagnetic radiation at the electron plasma frequency. This 
result can be used to offer an alternative interpretation to the horizontal 
strips (usually referred to as the narrowband emission lines in the literature) 
observed in some dynamical spectra. Quasilinear theory predictions: (i) 
electron free streaming and (ii) long relaxation time, in accord with the 
analytic expressions, are confirmed via direct, fully-kinetic simulation. 
• We considered interplay of Larmor-drift instability and plasma emission 
mechanism by studying a dense electron beam in the Larmor-drift unstable 
(inhomogencous) plasma. We established that in this case the physical sys- 
tem is driven by the effects of the beam while Larmor-drift is not dominant. 
We also found significant deviation from the quasilinear theory which man- 
ifests itself in (i) fast quasilinear relaxation, (ii) disintegration of the beam, 
and (iii) generation of significant electron return current and ion heating. 

We would like to close with the comments in relation to the prospects of 
comparison of the numerical simulations presented here with the observations 
(e.g. with the dynamical spectra - 2D radio emission intensity plots where fre- 
quency is on y-axis and time on x-axis) . Let us base our argument on the level of 
generated E y (EM component) in the considered three cases (Sections 2.1-2.3). 
In the Larmor-drift unstable case without the electron beam according to Figure 
1(c) E y attains values around 5 x 10 -7 . In the uniform (Larmor-drift stable case) 
with a low density beam (see. Figure 4(c)), E y attains only 10 -7 (five times less). 
In the Larmor-drift unstable case with dense beam E y attains 9 x 10 -4 (Figure 
6(c)). Based on this we conclude at this stage we cannot produce numerical 
(synthetic) dynamical spectrum where frequency of the radiation would drop 
in time as the beam moves along the decreasing density profile. This is for 
two reasons (i) when we consider small density beam rib/n e = 5 x 10~ 6 (when 
we are in the quasilinear regime and all known facts about plasma emission 
mechanism apply) the beam effect is too small {E y = 10~ 7 ) compared to the 
Larmor-drift unstable case without the electron beam (E y = 5x 10~ 7 ), thus 
there is no point presenting results of Larmor-drift unstable case with weak 
beam (because Larmor-drift instability overwhelms the effect of the weak beam), 
(ii) In the Larmor-drift unstable case with dense beam (nb/n e = 5 x 10 -2 ), we 
have E y = 9 x 10~ 4 , and we do not see decrease of the emission frequency with 
time because the beam disintegrates in quasilinear time of r = 20w~ e 1 (i.e. beam 
has not enough time to slide down the decreasing density profile) and hence we 
cannot expect to see plasma emission mechanism in action in its usual form. 
In summary, the model presented here cannot be directly applied to the type 
III busts and it is more relevant for the interpretation of the narrowband line 
emission observations. In oder to simulate the dynamical spectra of type III 
bursts in which the emission intensity rapidly drifts towards small frequencies in 
time, as the beam moves to a plasma with decreasing density (and hence w pe ), 
we would first need to suppress the Larmor-drift instability. In turn, this can be 
achieved by satisfying the condition set out in Equation (16). However, with the 
VALIS numerical code this is not possible because it only solves for (E x , E y , 0) 
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and (0, 0, B z ). We plan to study this issue further in a following publication (work 
in progress), using EPOCH 1.5D particle-in-cell code which allows to choose all 
background magnetic field components. Naturally, PIC method will suffer from 
the known shortcomings compared to the superior Vlasov-Maxwell approach. 
However, the benefit of ability of specifying all EM field components outweighs 
the downsides of the PIC approach. 
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